From 0a97edc2cc0a5b9a27c9302c59b5273c7179b4ab Mon Sep 17 00:00:00 2001 From: Drew Parsons Date: Wed, 14 Nov 2018 23:10:51 +0800 Subject: [PATCH] Import pygalmesh_0.2.6.orig.tar.gz [dgit import orig pygalmesh_0.2.6.orig.tar.gz] --- .circleci/config.yml | 21 ++ .flake8 | 5 + .gitignore | 11 + .pylintrc | 8 + .travis.yml | 34 ++ CMakeLists.txt | 5 + LICENSE | 21 ++ MANIFEST.in | 6 + Makefile | 32 ++ README.md | 321 ++++++++++++++++ pygalmesh/__about__.py | 11 + pygalmesh/__init__.py | 70 ++++ setup.py | 84 +++++ src/CMakeLists.txt | 31 ++ src/domain.hpp | 490 ++++++++++++++++++++++++ src/generate.cpp | 144 +++++++ src/generate.hpp | 33 ++ src/generate_from_off.cpp | 99 +++++ src/generate_from_off.hpp | 28 ++ src/generate_surface_mesh.cpp | 99 +++++ src/generate_surface_mesh.hpp | 23 ++ src/polygon2d.hpp | 308 +++++++++++++++ src/primitives.hpp | 453 ++++++++++++++++++++++ src/pybind11.cpp | 265 +++++++++++++ test/test_pygalmesh.py | 680 ++++++++++++++++++++++++++++++++++ 25 files changed, 3282 insertions(+) create mode 100644 .circleci/config.yml create mode 100644 .flake8 create mode 100644 .gitignore create mode 100644 .pylintrc create mode 100644 .travis.yml create mode 100644 CMakeLists.txt create mode 100644 LICENSE create mode 100644 MANIFEST.in create mode 100644 Makefile create mode 100644 README.md create mode 100644 pygalmesh/__about__.py create mode 100644 pygalmesh/__init__.py create mode 100644 setup.py create mode 100644 src/CMakeLists.txt create mode 100644 src/domain.hpp create mode 100644 src/generate.cpp create mode 100644 src/generate.hpp create mode 100644 src/generate_from_off.cpp create mode 100644 src/generate_from_off.hpp create mode 100644 src/generate_surface_mesh.cpp create mode 100644 src/generate_surface_mesh.hpp create mode 100644 src/polygon2d.hpp create mode 100644 src/primitives.hpp create mode 100644 src/pybind11.cpp create mode 100644 test/test_pygalmesh.py diff --git a/.circleci/config.yml b/.circleci/config.yml new file mode 100644 index 0000000..de4f763 --- /dev/null +++ b/.circleci/config.yml @@ -0,0 +1,21 @@ +version: 2 + +jobs: + + build: + working_directory: ~/pygalmesh + docker: + - image: ubuntu:18.04 + steps: + - run: apt-get update + - run: apt-get install -y libcgal-dev libeigen3-dev python3-pip clang++-6.0 + - run: pip3 install -U pytest pytest-cov black flake8 meshio pybind11 + - checkout + # install + # Use clang++ for its smaller memory footprint. + - run: CC="clang++-6.0" pip3 install . + # lint + - run: LC_ALL=C.UTF-8 black --check setup.py pygalmesh/ test/*.py + - run: flake8 setup.py pygalmesh/ test/*.py + # The actual test + - run: cd test/ && pytest diff --git a/.flake8 b/.flake8 new file mode 100644 index 0000000..c321e71 --- /dev/null +++ b/.flake8 @@ -0,0 +1,5 @@ +[flake8] +ignore = E203, E266, E501, W503 +max-line-length = 80 +max-complexity = 18 +select = B,C,E,F,W,T4,B9 diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..62bde14 --- /dev/null +++ b/.gitignore @@ -0,0 +1,11 @@ +*.off +*.vtu +*.mesh +.cache/ +build/ +dist/ +MANIFEST +README.rst +do-configure.sh +.pytest_cache/ +*.egg-info/ diff --git a/.pylintrc b/.pylintrc new file mode 100644 index 0000000..31f545d --- /dev/null +++ b/.pylintrc @@ -0,0 +1,8 @@ +[MESSAGES CONTROL] + +disable= + fixme, + invalid-name, + missing-docstring, + no-member, + too-few-public-methods diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..69ee8c6 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,34 @@ +dist: trusty + +language: python + +python: + - '2.7' + - '3.6' + +addons: + apt: + sources: + - sourceline: 'ppa:nschloe/cgal-backports' + - sourceline: 'ppa:nschloe/eigen-backports' + - sourceline: 'ppa:nschloe/swig-backports' + packages: + - libcgal-dev + - libcgal-qt5-11 # bug in cgal, should be installed automatically + - libeigen3-dev + - pandoc + - python-numpy + - swig + +# test dependencies +before_install: + - pip install meshio + +install: + - pip install pytest + - pip install . + +script: + - cd test + - curl https://raw.githubusercontent.com/libigl/libigl-unit-tests/master/data/elephant.off > elephant.off + - pytest diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..fb0710f --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,5 @@ +cmake_minimum_required(VERSION 3.0) + +project(pygalmesh CXX) + +add_subdirectory(src) diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..281c7b5 --- /dev/null +++ b/LICENSE @@ -0,0 +1,21 @@ +The MIT License (MIT) + +Copyright (c) 2016-2018 Nico Schlömer + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..c2fb0da --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1,6 @@ +include src/domain.hpp +include src/generate_from_off.hpp +include src/generate.hpp +include src/generate_surface_mesh.hpp +include src/polygon2d.hpp +include src/primitives.hpp diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..57e7287 --- /dev/null +++ b/Makefile @@ -0,0 +1,32 @@ +VERSION=$(shell python3 -c "import pygalmesh; print(pygalmesh.__version__)") + +default: + @echo "\"make publish\"?" + +tag: + @if [ "$(shell git rev-parse --abbrev-ref HEAD)" != "master" ]; then exit 1; fi + @echo "Tagging v$(VERSION)..." + git tag v$(VERSION) + git push --tags + +upload: setup.py + # Make sure we're on the master branch + @if [ "$(shell git rev-parse --abbrev-ref HEAD)" != "master" ]; then exit 1; fi + rm -rf dist/* + python3 setup.py sdist + twine upload dist/*.tar.gz + # HTTPError: 400 Client Error: Binary wheel 'pygalmesh-0.2.0-cp27-cp27mu-linux_x86_64.whl' has an unsupported platform tag 'linux_x86_64'. for url: https://upload.pypi.org/legacy/ + # python3 setup.py bdist_wheel --universal + # twine upload dist/*.whl + +publish: tag upload + +clean: + @find . | grep -E "(__pycache__|\.pyc|\.pyo$\)" | xargs rm -rf + +black: + black setup.py pygalmesh/ test/*.py + +lint: + black --check setup.py pygalmesh/ test/*.py + flake8 setup.py pygalmesh/ test/*.py diff --git a/README.md b/README.md new file mode 100644 index 0000000..fc45727 --- /dev/null +++ b/README.md @@ -0,0 +1,321 @@ +# pygalmesh + +A Python frontend to [CGAL](https://www.cgal.org/)'s [3D mesh generation +capabilities](https://doc.cgal.org/latest/Mesh_3/index.html). + +[![CircleCI](https://img.shields.io/circleci/project/github/nschloe/pygalmesh/master.svg)](https://circleci.com/gh/nschloe/pygalmesh/tree/master) +[![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/ambv/black) +[![PyPi Version](https://img.shields.io/pypi/v/pygalmesh.svg)](https://pypi.org/project/pygalmesh) +[![GitHub stars](https://img.shields.io/github/stars/nschloe/pygalmesh.svg?label=Stars&logo=github)](https://github.com/nschloe/pygalmesh) + +pygalmesh makes it easy to create high-quality 3D volume and surface meshes. + +### Background + +CGAL offers two different approaches for mesh generation: + +1. Meshes defined implicitly by level sets of functions. +2. Meshes defined by a set of bounding planes. + +pygalmesh provides a front-end to the first approach, which has the following +advantages and disadvantages: + +* All boundary points are guaranteed to be in the level set within any specified + residual. This results in smooth curved surfaces. +* Sharp intersections of subdomains (e.g., in unions or differences of sets) + need to be specified manually (via feature edges, see below), which can be + tedious. + +On the other hand, the bounding-plane approach (realized by +[mshr](https://bitbucket.org/fenics-project/mshr)), has the following +properties: + +* Smooth, curved domains are approximated by a set of bounding planes, + resulting in more of less visible edges. +* Intersections of domains can be computed automatically, so domain unions etc. + have sharp edges where they belong. + +Other Python mesh generators are [pygmsh](https://github.com/nschloe/pygmsh) (a +frontend to [gmsh](http://gmsh.info/)) and +[MeshPy](https://github.com/inducer/meshpy). +[meshzoo](https://github.com/nschloe/meshzoo) provides some basic canonical +meshes. + +### Examples + +#### A simple ball + + +```python +import pygalmesh + +s = pygalmesh.Ball([0, 0, 0], 1.0) +pygalmesh.generate_mesh(s, "out.mesh", cell_size=0.2) +``` +CGAL's mesh generator returns Medit-files, which can be processed by, e.g., +[meshio](https://github.com/nschloe/meshio). +```python +import meshio +vertices, cells, _, _, _ = meshio.read("out.mesh") +``` +The mesh generation comes with many more options, described +[here](https://doc.cgal.org/latest/Mesh_3/). Try, for example, +```python +pygalmesh.generate_mesh( + s, + "out.mesh", + cell_size=0.2, + edge_size=0.1, + odt=True, + lloyd=True, + verbose=False +) +``` + +#### Other primitive shapes + + +pygalmesh provides out-of-the-box support for balls, cuboids, ellipsoids, tori, +cones, cylinders, and tetrahedra. Try for example +```python +import pygalmesh + +s0 = pygalmesh.Tetrahedron( + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 1.0] +) +pygalmesh.generate_mesh(s0, "out.mesh", cell_size=0.1, edge_size=0.1) +``` + +#### Domain combinations + + +Supported are unions, intersections, and differences of all domains. As +mentioned above, however, the sharp intersections between two domains are not +automatically handled. Try for example +```python +import pygalmesh + +radius = 1.0 +displacement = 0.5 +s0 = pygalmesh.Ball([displacement, 0, 0], radius) +s1 = pygalmesh.Ball([-displacement, 0, 0], radius) +u = pygalmesh.Difference(s0, s1) +``` +To sharpen the intersection circle, add it as a feature edge polygon line, +e.g., +```python +a = numpy.sqrt(radius**2 - displacement**2) +edge_size = 0.15 +n = int(2*numpy.pi*a / edge_size) +circ = [ + [ + 0.0, + a * numpy.cos(i * 2*numpy.pi / n), + a * numpy.sin(i * 2*numpy.pi / n) + ] for i in range(n) + ] +circ.append(circ[0]) + +pygalmesh.generate_mesh( + u, + "out.mesh", + feature_edges=[circ], + cell_size=0.15, + edge_size=edge_size, + facet_angle=25, + facet_size=0.15, + cell_radius_edge_ratio=2.0 +) +``` +Note that the length of the polygon legs are kept in sync with the `edge_size` +of the mesh generation. This makes sure that it fits in nicely with the rest of +the mesh. + +#### Domain deformations + + +You can of course translate, rotate, scale, and stretch any domain. Try, for +example, +```python +import pygalmesh + +s = pygalmesh.Stretch( + pygalmesh.Ball([0, 0, 0], 1.0), + [1.0, 2.0, 0.0] +) + +pygalmesh.generate_mesh(s, "out.mesh", cell_size=0.1) +``` + +#### Extrusion of 2D polygons + + +pygalmesh lets you extrude any polygon into a 3D body. It even supports +rotation alongside! +```python +import pygalmesh + +p = pygalmesh.Polygon2D([[-0.5, -0.3], [0.5, -0.3], [0.0, 0.5]]) +edge_size = 0.1 +domain = pygalmesh.Extrude( + p, + [0.0, 0.0, 1.0], + 0.5 * 3.14159265359, + edge_size +) +pygalmesh.generate_mesh( + domain, + "out.mesh", + cell_size=0.1, + edge_size=edge_size, + verbose=False +) +``` +Feature edges are automatically preserved here, which is why an edge length +needs to be given to `pygalmesh.Extrude`. + +#### Rotation bodies + + +Polygons in the x-z-plane can also be rotated around the z-axis to yield a +rotation body. +```python +import pygalmesh + +p = pygalmesh.Polygon2D([[0.5, -0.3], [1.5, -0.3], [1.0, 0.5]]) +edge_size = 0.1 +domain = pygalmesh.RingExtrude(p, edge_size) +pygalmesh.generate_mesh( + domain, + "out.mesh", + cell_size=0.1, + edge_size=edge_size, + verbose=False +) +``` + +#### Your own custom level set function + + +If all of the variety is not enough for you, you can define your own custom +level set function. You simply need to subclass `pygalmesh.DomainBase` and +specify a function, e.g., +```python +import pygalmesh +class Heart(pygalmesh.DomainBase): + def __init__(self): + super(Heart, self).__init__() + return + + def eval(self, x): + return (x[0]**2 + 9.0/4.0 * x[1]**2 + x[2]**2 - 1)**3 \ + - x[0]**2 * x[2]**3 - 9.0/80.0 * x[1]**2 * x[2]**3 + + def get_bounding_sphere_squared_radius(self): + return 10.0 + +d = Heart() +pygalmesh.generate_mesh(d, "out.mesh", cell_size=0.1) +``` +Note that you need to specify the square of a bounding sphere radius, used as +an input to CGAL's mesh generator. + +#### Surface meshes + +If you're only after the surface of a body, pygalmesh has +`generate_surface_mesh` for you. It offers fewer options (obviously, +`cell_size` is gone), but otherwise works the same way: +```python +import pygalmesh + +s = pygalmesh.Ball([0, 0, 0], 1.0) +pygalmesh.generate_surface_mesh( + s, + "out.off", + angle_bound=30, + radius_bound=0.1, + distance_bound=0.1 +) +``` +The output format is +[OFF](http://segeval.cs.princeton.edu/public/off_format.html) which again is +handled by [meshio](https://github.com/nschloe/meshio). + +Refer to [CGAL's +documention](https://doc.cgal.org/latest/Surface_mesher/index.html) for the +options. + +#### Meshes from OFF files + + +If you have an OFF file at hand (like +[elephant.off](https://raw.githubusercontent.com/CGAL/cgal-swig-bindings/master/examples/data/elephant.off) +or [these](https://github.com/CGAL/cgal/tree/master/Surface_mesher/demo/Surface_mesher/inputs)), +pygalmesh generates the mesh via +```python +import pygalmesh + +pygalmesh.generate_from_off( + "elephant.off", + "out.mesh", + facet_angle=25.0, + facet_size=0.15, + facet_distance=0.008, + cell_radius_edge_ratio=3.0, + verbose=False +) +``` + +### Installation + +For installation, pygalmesh needs [CGAL](https://www.cgal.org/) and +[Eigen](http://eigen.tuxfamily.org/index.php?title=Main_Page) installed on your +system. They are typically available on your Linux distribution, e.g., on +Ubuntu +``` +sudo apt install libcgal-dev libeigen3-dev +``` +After that, pygalmesh can be [installed from the Python Package +Index](https://pypi.org/project/pygalmesh/), so with +``` +pip install -U pygalmesh +``` +you can install/upgrade. + +[meshio](https://github.com/nschloe/meshio) (`sudo -H pip install meshio`) +can be helpful in processing the meshes. + +#### Manual installation + +For manual installation (if you're a developer or just really keen on getting +the bleeding edge version of pygalmesh), there are two possibilities: + + * Get the sources, type `sudo python setup.py install`. This does the trick + most the time. + * As a fallback, there's a CMake-based installation. Simply go `cmake + /path/to/sources/` and `make`. + +### Testing + +To run the pygalmesh unit tests, check out this repository and type +``` +pytest +``` + +### Distribution + +To create a new release + +1. bump the `__version__` number (in `setup.py` _and_ `src/pygalmesh.i`) + +2. publish to PyPi and GitHub: + ``` + make publish + ``` + +### License + +pygalmesh is published under the [MIT license](https://en.wikipedia.org/wiki/MIT_License). diff --git a/pygalmesh/__about__.py b/pygalmesh/__about__.py new file mode 100644 index 0000000..5bf2aa0 --- /dev/null +++ b/pygalmesh/__about__.py @@ -0,0 +1,11 @@ +# -*- coding: utf-8 -*- +# + +__author__ = u"Nico Schlömer" +__author_email__ = "nico.schloemer@gmail.com" +__copyright__ = u"Copyright (c) 2016-2018, {} <{}>".format(__author__, __author_email__) +__license__ = "License :: OSI Approved :: MIT License" +__version__ = "0.2.6" +__maintainer__ = u"Nico Schlömer" +__status__ = "Development Status :: 4 - Beta" +__url__ = "https://github.com/nschloe/pygalmesh" diff --git a/pygalmesh/__init__.py b/pygalmesh/__init__.py new file mode 100644 index 0000000..dc00c34 --- /dev/null +++ b/pygalmesh/__init__.py @@ -0,0 +1,70 @@ +# -*- coding: utf-8 -*- +# +# https://github.com/pybind/pybind11/issues/1004 +from _pygalmesh import ( + DomainBase, + Translate, + Rotate, + Scale, + Stretch, + Intersection, + Union, + Difference, + Extrude, + Ball, + Cuboid, + Ellipsoid, + Tetrahedron, + Cone, + Cylinder, + Torus, + HalfSpace, + Polygon2D, + RingExtrude, + generate_mesh, + generate_surface_mesh, + generate_from_off, +) + +from .__about__ import ( + __author__, + __author_email__, + __copyright__, + __license__, + __version__, + __maintainer__, + __status__, +) + +__all__ = [ + "__author__", + "__author_email__", + "__copyright__", + "__license__", + "__version__", + "__maintainer__", + "__status__", + # + "DomainBase", + "Translate", + "Rotate", + "Scale", + "Stretch", + "Intersection", + "Union", + "Difference", + "Extrude", + "Ball", + "Cuboid", + "Ellipsoid", + "Tetrahedron", + "Cone", + "Cylinder", + "Torus", + "HalfSpace", + "Polygon2D", + "RingExtrude", + "generate_mesh", + "generate_surface_mesh", + "generate_from_off", +] diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..16663a8 --- /dev/null +++ b/setup.py @@ -0,0 +1,84 @@ +# -*- coding: utf-8 -*- +# +import os + +import codecs +from setuptools import setup, Extension, find_packages + + +# https://packaging.python.org/single_source_version/ +base_dir = os.path.abspath(os.path.dirname(__file__)) +about = {} +with open(os.path.join(base_dir, "pygalmesh", "__about__.py"), "rb") as handle: + # pylint: disable=exec-used + exec(handle.read(), about) + + +class get_pybind_include(object): + """Helper class to determine the pybind11 include path + The purpose of this class is to postpone importing pybind11 + until it is actually installed, so that the ``get_include()`` + method can be invoked. + """ + + def __init__(self, user=False): + self.user = user + + def __str__(self): + import pybind11 + + return pybind11.get_include(self.user) + + +def read(fname): + return codecs.open(os.path.join(base_dir, fname), encoding="utf-8").read() + + +ext_modules = [ + Extension( + "_pygalmesh", + [ + "src/generate.cpp", + "src/generate_from_off.cpp", + "src/generate_surface_mesh.cpp", + "src/pybind11.cpp", + ], + language="c++", + include_dirs=[ + "/usr/include/eigen3/", + # Path to pybind11 headers + get_pybind_include(), + get_pybind_include(user=True), + ], + libraries=["CGAL", "gmp", "mpfr"], + # extra_compile_args=['-std=c++11'] + ) +] + +setup( + name="pygalmesh", + packages=find_packages(), + # cmdclass={'build_ext': BuildExt}, + ext_modules=ext_modules, + version=about["__version__"], + url=about["__url__"], + author=about["__author__"], + author_email=about["__author_email__"], + install_requires=["numpy", "pybind11 >= 2.2", "pipdate"], + description="Python frontend to CGAL's 3D mesh generation capabilities", + long_description=read("README.md"), + long_description_content_type="text/markdown", + license=about["__license__"], + classifiers=[ + about["__status__"], + about["__license__"], + "Operating System :: OS Independent", + "Programming Language :: Python", + "Programming Language :: Python :: 2", + "Programming Language :: Python :: 3", + "Topic :: Scientific/Engineering", + "Topic :: Scientific/Engineering :: Mathematics", + "Topic :: Scientific/Engineering :: Physics", + "Topic :: Scientific/Engineering :: Visualization", + ], +) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt new file mode 100644 index 0000000..a6df64b --- /dev/null +++ b/src/CMakeLists.txt @@ -0,0 +1,31 @@ +FIND_PACKAGE(pybind11 REQUIRED) +# include_directories(${PYBIND11_INCLUDE_DIR}) + +FIND_PACKAGE(Eigen3 REQUIRED) +include_directories(${EIGEN3_INCLUDE_DIR}) + +FILE(GLOB pygalmesh_SRCS "${CMAKE_CURRENT_SOURCE_DIR}/*.cpp") +# FILE(GLOB pygalmesh_HEADERS "${CMAKE_CURRENT_SOURCE_DIR}/*.hpp") + +pybind11_add_module(pygalmesh ${pygalmesh_SRCS}) + +# Call CGAL as late as possible. It messes around with some +# compiler variables. See bugs +# and +# . +FIND_PACKAGE(CGAL REQUIRED) +include(${CGAL_USE_FILE}) + +# ADD_LIBRARY(pygalmesh ${pygalmesh_SRCS}) +target_link_libraries(pygalmesh PRIVATE ${CGAL_LIBRARIES}) + +# execute_process( +# COMMAND python -c "from distutils.sysconfig import get_python_lib; print get_python_lib()" +# OUTPUT_VARIABLE PYTHON_SITE_PACKAGES +# OUTPUT_STRIP_TRAILING_WHITESPACE +# ) +# install(TARGETS _pygalmesh DESTINATION ${PYTHON_SITE_PACKAGES}) +# install( +# FILES ${CMAKE_BINARY_DIR}/src/pygalmesh.py +# DESTINATION ${PYTHON_SITE_PACKAGES} +# ) diff --git a/src/domain.hpp b/src/domain.hpp new file mode 100644 index 0000000..ff8cd9c --- /dev/null +++ b/src/domain.hpp @@ -0,0 +1,490 @@ +#ifndef DOMAIN_HPP +#define DOMAIN_HPP + +#include +#include +#include +#include +#include + +namespace pygalmesh { + +class DomainBase +{ + public: + + virtual ~DomainBase() = default; + + virtual + double + eval(const std::array & x) const = 0; + + virtual + double + get_bounding_sphere_squared_radius() const = 0; + + virtual + std::vector>> + get_features() const + { + return {}; + }; +}; + +class Translate: public pygalmesh::DomainBase +{ + public: + Translate( + const std::shared_ptr & domain, + const std::array & direction + ): + domain_(domain), + direction_(Eigen::Vector3d(direction.data())), + translated_features_(translate_features(domain->get_features(), direction_)) + { + } + + virtual ~Translate() = default; + + std::vector>> + translate_features( + const std::vector>> & features, + const Eigen::Vector3d & direction + ) const + { + std::vector>> translated_features; + for (const auto & feature: features) { + std::vector> translated_feature; + for (const auto & point: feature) { + const std::array translated_point = { + point[0] + direction[0], + point[1] + direction[1], + point[2] + direction[2] + }; + translated_feature.push_back(translated_point); + } + translated_features.push_back(translated_feature); + } + return translated_features; + } + + virtual + double + eval(const std::array & x) const + { + const std::array d = { + x[0] - direction_[0], + x[1] - direction_[1], + x[2] - direction_[2] + }; + return domain_->eval(d); + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + const double radius = sqrt(domain_->get_bounding_sphere_squared_radius()); + const double dir_norm = direction_.norm(); + return (radius + dir_norm)*(radius + dir_norm); + } + + virtual + std::vector>> + get_features() const + { + return translated_features_; + }; + + private: + const std::shared_ptr domain_; + const Eigen::Vector3d direction_; + const std::vector>> translated_features_; +}; + +class Rotate: public pygalmesh::DomainBase +{ + public: + Rotate( + const std::shared_ptr & domain, + const std::array & axis, + const double angle + ): + domain_(domain), + normalized_axis_(Eigen::Vector3d(axis.data()).normalized()), + sinAngle_(sin(angle)), + cosAngle_(cos(angle)), + rotated_features_(rotate_features(domain_->get_features())) + { + } + + virtual ~Rotate() = default; + + Eigen::Vector3d + rotate( + const Eigen::Vector3d & vec, + const Eigen::Vector3d & axis, + const double sinAngle, + const double cosAngle + ) const + { + // Rotate a vector `v` by the angle `theta` in the plane perpendicular + // to the axis given by `u`. + // Refer to + // http://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle + // + // cos(theta) * I * v + // + sin(theta) u\cross v + // + (1-cos(theta)) (u*u^T) * v + return cosAngle * vec + + sinAngle * axis.cross(vec) + + (1.0-cosAngle) * axis.dot(vec) * axis; + } + + virtual + double + eval(const std::array & x) const + { + // rotate with negative angle + const auto p2 = rotate( + Eigen::Vector3d(x.data()), + normalized_axis_, + -sinAngle_, + cosAngle_ + ); + return domain_->eval({p2[0], p2[1], p2[2]}); + } + + std::vector>> + rotate_features( + const std::vector>> & features + ) const + { + std::vector>> rotated_features; + for (const auto & feature: features) { + std::vector> rotated_feature; + for (const auto & point: feature) { + const auto p2 = rotate( + Eigen::Vector3d(point.data()), + normalized_axis_, + sinAngle_, + cosAngle_ + ); + rotated_feature.push_back({p2[0], p2[1], p2[2]}); + } + rotated_features.push_back(rotated_feature); + } + return rotated_features; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + return domain_->get_bounding_sphere_squared_radius(); + } + + virtual + std::vector>> + get_features() const + { + return rotated_features_; + }; + + private: + const std::shared_ptr domain_; + const Eigen::Vector3d normalized_axis_; + const double sinAngle_; + const double cosAngle_; + const std::vector>> rotated_features_; +}; + +class Scale: public pygalmesh::DomainBase +{ + public: + Scale( + std::shared_ptr & domain, + const double alpha + ): + domain_(domain), + alpha_(alpha), + scaled_features_(scale_features(domain_->get_features())) + { + assert(alpha_ > 0.0); + } + + virtual ~Scale() = default; + + virtual + double + eval(const std::array & x) const + { + return domain_->eval({x[0]/alpha_, x[1]/alpha_, x[2]/alpha_}); + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + return alpha_*alpha_ * domain_->get_bounding_sphere_squared_radius(); + } + + std::vector>> + scale_features( + const std::vector>> & features + ) const + { + std::vector>> scaled_features; + for (const auto & feature: features) { + std::vector> scaled_feature; + for (const auto & point: feature) { + scaled_feature.push_back({ + alpha_ * point[0], + alpha_ * point[1], + alpha_ * point[2] + }); + } + scaled_features.push_back(scaled_feature); + } + return scaled_features; + } + + virtual + std::vector>> + get_features() const + { + return scaled_features_; + }; + + private: + std::shared_ptr domain_; + const double alpha_; + const std::vector>> scaled_features_; +}; + +class Stretch: public pygalmesh::DomainBase +{ + public: + Stretch( + std::shared_ptr & domain, + const std::array & direction + ): + domain_(domain), + normalized_direction_(Eigen::Vector3d(direction.data()).normalized()), + alpha_(Eigen::Vector3d(direction.data()).norm()), + stretched_features_(stretch_features(domain_->get_features())) + { + assert(alpha_ > 0.0); + } + + virtual ~Stretch() = default; + + virtual + double + eval(const std::array & x) const + { + const Eigen::Vector3d v(x.data()); + const double beta = normalized_direction_.dot(v); + // scale the component of normalized_direction_ by 1/alpha_ + const auto v2 = beta/alpha_ * normalized_direction_ + + (v - beta * normalized_direction_); + return domain_->eval({v2[0], v2[1], v2[2]}); + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + return alpha_*alpha_ * domain_->get_bounding_sphere_squared_radius(); + } + + std::vector>> + stretch_features( + const std::vector>> & features + ) const + { + std::vector>> stretched_features; + for (const auto & feature: features) { + std::vector> stretched_feature; + for (const auto & point: feature) { + // scale the component of normalized_direction_ by alpha_ + const Eigen::Vector3d v(point.data()); + const double beta = normalized_direction_.dot(v); + const auto v2 = beta * alpha_ * normalized_direction_ + + (v - beta * normalized_direction_); + stretched_feature.push_back({v2[0], v2[1], v2[2]}); + } + stretched_features.push_back(stretched_feature); + } + return stretched_features; + } + + virtual + std::vector>> + get_features() const + { + return stretched_features_; + }; + + private: + std::shared_ptr domain_; + const Eigen::Vector3d normalized_direction_; + const double alpha_; + const std::vector>> stretched_features_; +}; + +class Intersection: public pygalmesh::DomainBase +{ + public: + explicit Intersection( + std::vector> & domains + ): + domains_(domains) + { + } + + virtual ~Intersection() = default; + + virtual + double + eval(const std::array & x) const + { + // TODO find a differentiable expression + double maxval = std::numeric_limits::lowest(); + for (const auto & domain: domains_) { + maxval = std::max(maxval, domain->eval(x)); + } + return maxval; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + double min = std::numeric_limits::max(); + for (const auto & domain: domains_) { + min = std::min(min, domain->get_bounding_sphere_squared_radius()); + } + return min; + } + + virtual + std::vector>> + get_features() const + { + std::vector>> features; + for (const auto & domain: domains_) { + const auto f = domain->get_features(); + features.insert(std::end(features), std::begin(f), std::end(f)); + } + return features; + }; + + private: + std::vector> domains_; +}; + +class Union: public pygalmesh::DomainBase +{ + public: + explicit Union( + std::vector> & domains + ): + domains_(domains) + { + } + + virtual ~Union() = default; + + virtual + double + eval(const std::array & x) const + { + // TODO find a differentiable expression + double minval = std::numeric_limits::max(); + for (const auto & domain: domains_) { + minval = std::min(minval, domain->eval(x)); + } + return minval; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + double max = 0.0; + for (const auto & domain: domains_) { + max = std::max(max, domain->get_bounding_sphere_squared_radius()); + } + return max; + } + + virtual + std::vector>> + get_features() const + { + std::vector>> features; + for (const auto & domain: domains_) { + const auto f = domain->get_features(); + features.insert(std::end(features), std::begin(f), std::end(f)); + } + return features; + }; + + private: + std::vector> domains_; +}; + +class Difference: public pygalmesh::DomainBase +{ + public: + Difference( + std::shared_ptr & domain0, + std::shared_ptr & domain1 + ): + domain0_(domain0), + domain1_(domain1) + { + } + + virtual ~Difference() = default; + + virtual + double + eval(const std::array & x) const + { + // TODO find a continuous (perhaps even differentiable) expression + const double val0 = domain0_->eval(x); + const double val1 = domain1_->eval(x); + return (val0 < 0.0 && val1 >= 0.0) ? val0 : std::max(val0, -val1); + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + return domain0_->get_bounding_sphere_squared_radius(); + } + + virtual + std::vector>> + get_features() const + { + std::vector>> features; + + const auto f0 = domain0_->get_features(); + features.insert(std::end(features), std::begin(f0), std::end(f0)); + + const auto f1 = domain1_->get_features(); + features.insert(std::end(features), std::begin(f1), std::end(f1)); + + return features; + }; + + private: + std::shared_ptr domain0_; + std::shared_ptr domain1_; +}; + +} // namespace pygalmesh +#endif // DOMAIN_HPP diff --git a/src/generate.cpp b/src/generate.cpp new file mode 100644 index 0000000..0e2e45f --- /dev/null +++ b/src/generate.cpp @@ -0,0 +1,144 @@ +#define CGAL_MESH_3_VERBOSE 1 + +#include "generate.hpp" + +#include + +#include +#include +#include + +#include +#include +#include + +namespace pygalmesh { + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; + +// Wrapper for DomainBase for translating to K::Point_3/FT. +class CgalDomainWrapper +{ + public: + explicit CgalDomainWrapper(const std::shared_ptr & domain): + domain_(domain) + { + } + + virtual ~CgalDomainWrapper() = default; + + virtual + K::FT + operator()(K::Point_3 p) const + { + return domain_->eval({p.x(), p.y(), p.z()}); + } + + private: + const std::shared_ptr domain_; +}; + +typedef CGAL::Mesh_domain_with_polyline_features_3< + CGAL::Implicit_mesh_domain_3 + > + Mesh_domain; + +// Triangulation +typedef CGAL::Mesh_triangulation_3::type Tr; +typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3; + +// Mesh Criteria +typedef CGAL::Mesh_criteria_3 Mesh_criteria; +typedef Mesh_criteria::Facet_criteria Facet_criteria; +typedef Mesh_criteria::Cell_criteria Cell_criteria; + +// translate vector> to list> +std::list> +translate_feature_edges( + const std::vector>> & feature_edges + ) +{ + std::list> polylines; + for (const auto & feature_edge: feature_edges) { + std::vector polyline; + for (const auto & point: feature_edge) { + polyline.push_back(K::Point_3(point[0], point[1], point[2])); + } + polylines.push_back(polyline); + } + return polylines; +} + +void +generate_mesh( + const std::shared_ptr & domain, + const std::string & outfile, + const std::vector>> & feature_edges, + const double bounding_sphere_radius, + const double boundary_precision, + const bool lloyd, + const bool odt, + const bool perturb, + const bool exude, + const double edge_size, + const double facet_angle, + const double facet_size, + const double facet_distance, + const double cell_radius_edge_ratio, + const double cell_size, + const bool verbose + ) +{ + const double bounding_sphere_radius2 = bounding_sphere_radius > 0 ? + bounding_sphere_radius*bounding_sphere_radius : + // some wiggle room + 1.01 * domain->get_bounding_sphere_squared_radius(); + + const auto d = CgalDomainWrapper(domain); + Mesh_domain cgal_domain( + d, + K::Sphere_3(CGAL::ORIGIN, bounding_sphere_radius2), + boundary_precision + ); + + const auto native_features = translate_feature_edges(domain->get_features()); + cgal_domain.add_features(native_features.begin(), native_features.end()); + + const auto polylines = translate_feature_edges(feature_edges); + cgal_domain.add_features(polylines.begin(), polylines.end()); + + Mesh_criteria criteria( + CGAL::parameters::edge_size=edge_size, + CGAL::parameters::facet_angle=facet_angle, + CGAL::parameters::facet_size=facet_size, + CGAL::parameters::facet_distance=facet_distance, + CGAL::parameters::cell_radius_edge_ratio=cell_radius_edge_ratio, + CGAL::parameters::cell_size=cell_size + ); + + // Mesh generation + if (!verbose) { + // suppress output + std::cerr.setstate(std::ios_base::failbit); + } + C3t3 c3t3 = CGAL::make_mesh_3( + cgal_domain, + criteria, + lloyd ? CGAL::parameters::lloyd() : CGAL::parameters::no_lloyd(), + odt ? CGAL::parameters::odt() : CGAL::parameters::no_odt(), + perturb ? CGAL::parameters::perturb() : CGAL::parameters::no_perturb(), + exude ? CGAL::parameters::exude() : CGAL::parameters::no_exude() + ); + if (!verbose) { + std::cerr.clear(); + } + + // Output + std::ofstream medit_file(outfile); + c3t3.output_to_medit(medit_file); + medit_file.close(); + + return; +} + +} // namespace pygalmesh diff --git a/src/generate.hpp b/src/generate.hpp new file mode 100644 index 0000000..f49fa31 --- /dev/null +++ b/src/generate.hpp @@ -0,0 +1,33 @@ +#ifndef GENERATE_HPP +#define GENERATE_HPP + +#include "domain.hpp" + +#include +#include +#include + +namespace pygalmesh { + +void generate_mesh( + const std::shared_ptr & domain, + const std::string & outfile, + const std::vector>> & feature_edges = {}, + const double bounding_sphere_radius = 0.0, + const double boundary_precision = 1.0e-4, + const bool lloyd = false, + const bool odt = false, + const bool perturb = true, + const bool exude = true, + const double edge_size = 0.0, // std::numeric_limits::max(), + const double facet_angle = 0.0, + const double facet_size = 0.0, + const double facet_distance = 0.0, + const double cell_radius_edge_ratio = 0.0, + const double cell_size = 0.0, + const bool verbose = true + ); + +} // namespace pygalmesh + +#endif // GENERATE_HPP diff --git a/src/generate_from_off.cpp b/src/generate_from_off.cpp new file mode 100644 index 0000000..83ab7e9 --- /dev/null +++ b/src/generate_from_off.cpp @@ -0,0 +1,99 @@ +#include "generate_from_off.hpp" + +#include +#include +#include +#include +#include +#include +#include + +// IO +#include + +namespace pygalmesh { + +// Domain +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +typedef CGAL::Polyhedron_3 Polyhedron; +typedef CGAL::Polyhedral_mesh_domain_3 Mesh_domain; + +// Triangulation +typedef CGAL::Mesh_triangulation_3::type Tr; + +typedef CGAL::Mesh_complex_3_in_triangulation_3 C3t3; + +// Criteria +typedef CGAL::Mesh_criteria_3 Mesh_criteria; + +// To avoid verbose function and named parameters call +using namespace CGAL::parameters; + +void +generate_from_off( + const std::string & infile, + const std::string & outfile, + const bool lloyd, + const bool odt, + const bool perturb, + const bool exude, + const double edge_size, + const double facet_angle, + const double facet_size, + const double facet_distance, + const double cell_radius_edge_ratio, + const double cell_size, + const bool verbose + ) +{ + // Create input polyhedron + Polyhedron polyhedron; + std::ifstream input(infile); + input >> polyhedron; + if (!input.good()) { + std::stringstream msg; + msg << "Cannot read input file \"" << infile << "\"" << std::endl; + throw std::runtime_error(msg.str()); + } + input.close(); + + // Create domain + Mesh_domain cgal_domain(polyhedron); + + // Mesh criteria + Mesh_criteria criteria( + CGAL::parameters::edge_size=edge_size, + CGAL::parameters::facet_angle=facet_angle, + CGAL::parameters::facet_size=facet_size, + CGAL::parameters::facet_distance=facet_distance, + CGAL::parameters::cell_radius_edge_ratio=cell_radius_edge_ratio, + CGAL::parameters::cell_size=cell_size + ); + + // Mesh generation + if (!verbose) { + // suppress output + std::cerr.setstate(std::ios_base::failbit); + } + // TODO make_mesh_3 segfaults on travis. hm. + C3t3 c3t3 = CGAL::make_mesh_3( + cgal_domain, + criteria, + lloyd ? CGAL::parameters::lloyd() : CGAL::parameters::no_lloyd(), + odt ? CGAL::parameters::odt() : CGAL::parameters::no_odt(), + perturb ? CGAL::parameters::perturb() : CGAL::parameters::no_perturb(), + exude ? CGAL::parameters::exude() : CGAL::parameters::no_exude() + ); + if (!verbose) { + std::cerr.clear(); + } + + // Output + std::ofstream medit_file(outfile); + c3t3.output_to_medit(medit_file); + medit_file.close(); + + return; +} + +} // namespace pygalmesh diff --git a/src/generate_from_off.hpp b/src/generate_from_off.hpp new file mode 100644 index 0000000..c2aba9c --- /dev/null +++ b/src/generate_from_off.hpp @@ -0,0 +1,28 @@ +#ifndef GENERATE_FROM_OFF_HPP +#define GENERATE_FROM_OFF_HPP + +#include +#include + +namespace pygalmesh { + +void +generate_from_off( + const std::string & infile, + const std::string & outfile, + const bool lloyd = false, + const bool odt = false, + const bool perturb = true, + const bool exude = true, + const double edge_size = 0.0, // std::numeric_limits::max(), + const double facet_angle = 0.0, + const double facet_size = 0.0, + const double facet_distance = 0.0, + const double cell_radius_edge_ratio = 0.0, + const double cell_size = 0.0, + const bool verbose = true + ); + +} // namespace pygalmesh + +#endif // GENERATE_FROM_OFF_HPP diff --git a/src/generate_surface_mesh.cpp b/src/generate_surface_mesh.cpp new file mode 100644 index 0000000..8a38b02 --- /dev/null +++ b/src/generate_surface_mesh.cpp @@ -0,0 +1,99 @@ +#define CGAL_SURFACE_MESHER_VERBOSE 1 + +#include "generate_surface_mesh.hpp" + +#include +#include +#include +#include +#include +#include + +namespace pygalmesh { + +// default triangulation for Surface_mesher +typedef CGAL::Surface_mesh_default_triangulation_3 Tr; +// c2t3 +typedef CGAL::Complex_2_in_triangulation_3 C2t3; +typedef Tr::Geom_traits GT; + +// Wrapper for DomainBase for translating to GT. +class CgalDomainWrapper +{ + public: + explicit CgalDomainWrapper(const std::shared_ptr & domain): + domain_(domain) + { + } + + virtual ~CgalDomainWrapper() = default; + + virtual + GT::FT + operator()(GT::Point_3 p) const + { + return domain_->eval({p.x(), p.y(), p.z()}); + } + + private: + const std::shared_ptr domain_; +}; + +typedef CGAL::Implicit_surface_3 Surface_3; + +void +generate_surface_mesh( + const std::shared_ptr & domain, + const std::string & outfile, + const double bounding_sphere_radius, + const double angle_bound, + const double radius_bound, + const double distance_bound, + const bool verbose + ) +{ + const double bounding_sphere_radius2 = bounding_sphere_radius > 0 ? + bounding_sphere_radius*bounding_sphere_radius : + // add a little wiggle room + 1.01 * domain->get_bounding_sphere_squared_radius(); + + Tr tr; // 3D-Delaunay triangulation + C2t3 c2t3 (tr); // 2D-complex in 3D-Delaunay triangulation + + const auto d = CgalDomainWrapper(domain); + Surface_3 surface( + d, + GT::Sphere_3(CGAL::ORIGIN, bounding_sphere_radius2) + ); + + CGAL::Surface_mesh_default_criteria_3 criteria( + angle_bound, + radius_bound, + distance_bound + ); + + if (!verbose) { + // suppress output + std::cout.setstate(std::ios_base::failbit); + std::cerr.setstate(std::ios_base::failbit); + } + CGAL::make_surface_mesh( + c2t3, + surface, + criteria, + CGAL::Non_manifold_tag() + ); + if (!verbose) { + std::cout.clear(); + std::cerr.clear(); + } + + // Output + std::ofstream off_file(outfile); + CGAL::output_surface_facets_to_off(off_file, c2t3); + off_file.close(); + + return; +} + +} // namespace pygalmesh diff --git a/src/generate_surface_mesh.hpp b/src/generate_surface_mesh.hpp new file mode 100644 index 0000000..8a63043 --- /dev/null +++ b/src/generate_surface_mesh.hpp @@ -0,0 +1,23 @@ +#ifndef GENERATE_SURFACE_MESH_HPP +#define GENERATE_SURFACE_MESH_HPP + +#include "domain.hpp" + +#include +#include + +namespace pygalmesh { + +void generate_surface_mesh( + const std::shared_ptr & domain, + const std::string & outfile, + const double bounding_sphere_radius = 0.0, + const double angle_bound = 0.0, + const double radius_bound = 0.0, + const double distance_bound = 0.0, + const bool verbose = true + ); + +} // namespace pygalmesh + +#endif // GENERATE_SURFACE_MESH_HPP diff --git a/src/polygon2d.hpp b/src/polygon2d.hpp new file mode 100644 index 0000000..6a5e231 --- /dev/null +++ b/src/polygon2d.hpp @@ -0,0 +1,308 @@ +#ifndef POLYGON2D_HPP +#define POLYGON2D_HPP + +#include "domain.hpp" + +#include +#include +#include +#include +#include + +typedef CGAL::Exact_predicates_inexact_constructions_kernel K; + +namespace pygalmesh { + +class Polygon2D { + public: + explicit Polygon2D(const std::vector> & _points): + points(vector_to_cgal_points(_points)) + { + } + + virtual ~Polygon2D() = default; + + std::vector + vector_to_cgal_points(const std::vector> & _points) const + { + std::vector points2(_points.size()); + for (size_t i = 0; i < _points.size(); i++) { + assert(_points[i].size() == 2); + points2[i] = K::Point_2(_points[i][0], _points[i][1]); + } + return points2; + } + + bool + is_inside(const std::array & point) + { + K::Point_2 pt(point[0], point[1]); + switch(CGAL::bounded_side_2(this->points.begin(), this->points.end(), pt, K())) { + case CGAL::ON_BOUNDED_SIDE: + return true; + case CGAL::ON_BOUNDARY: + return true; + case CGAL::ON_UNBOUNDED_SIDE: + return false; + default: + return false; + } + return false; + } + + public: + const std::vector points; +}; + + +class Extrude: public pygalmesh::DomainBase { + public: + Extrude( + const std::shared_ptr & poly, + const std::array & direction, + const double alpha = 0.0, + const double edge_size = 0.0 + ): + poly_(poly), + direction_(direction), + alpha_(alpha), + edge_size_(edge_size) + { + } + + virtual ~Extrude() = default; + + virtual + double + eval(const std::array & x) const + { + if (x[2] < 0.0 || x[2] > direction_[2]) { + return 1.0; + } + + const double beta = x[2] / direction_[2]; + + std::array x2 = { + x[0] - beta * direction_[0], + x[1] - beta * direction_[1] + }; + + if (alpha_ != 0.0) { + std::array x3; + // turn by -beta*alpha + const double sinAlpha = sin(beta*alpha_); + const double cosAlpha = cos(beta*alpha_); + x3[0] = cosAlpha * x2[0] + sinAlpha * x2[1]; + x3[1] = -sinAlpha * x2[0] + cosAlpha * x2[1]; + x2 = x3; + } + + return poly_->is_inside(x2) ? -1.0 : 1.0; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + double max = 0.0; + for (const auto & pt: poly_->points) { + // bottom polygon + const double nrm0 = pt.x()*pt.x() + pt.y()*pt.y(); + if (nrm0 > max) { + max = nrm0; + } + + // TODO rotation + + // top polygon + const double x = pt.x() + direction_[0]; + const double y = pt.y() + direction_[1]; + const double z = direction_[2]; + const double nrm1 = x*x + y*y + z*z; + if (nrm1 > max) { + max = nrm1; + } + } + return max; + } + + virtual + std::vector>> + get_features() const + { + std::vector>> features = {}; + + size_t n; + + // bottom polygon + n = poly_->points.size(); + for (size_t i=0; i < n-1; i++) { + features.push_back({ + {poly_->points[i].x(), poly_->points[i].y(), 0.0}, + {poly_->points[i+1].x(), poly_->points[i+1].y(), 0.0} + }); + } + features.push_back({ + {poly_->points[n-1].x(), poly_->points[n-1].y(), 0.0}, + {poly_->points[0].x(), poly_->points[0].y(), 0.0} + }); + + // top polygon, R*x + d + n = poly_->points.size(); + const double sinAlpha = sin(alpha_); + const double cosAlpha = cos(alpha_); + for (size_t i=0; i < n-1; i++) { + features.push_back({ + { + cosAlpha * poly_->points[i].x() - sinAlpha * poly_->points[i].y() + direction_[0], + sinAlpha * poly_->points[i].x() + cosAlpha * poly_->points[i].y() + direction_[1], + direction_[2] + }, + { + cosAlpha * poly_->points[i+1].x() - sinAlpha * poly_->points[i+1].y() + direction_[0], + sinAlpha * poly_->points[i+1].x() + cosAlpha * poly_->points[i+1].y() + direction_[1], + direction_[2] + } + }); + } + features.push_back({ + { + cosAlpha * poly_->points[n-1].x() - sinAlpha * poly_->points[n-1].y() + direction_[0], + sinAlpha * poly_->points[n-1].x() + cosAlpha * poly_->points[n-1].y() + direction_[1], + direction_[2] + }, + { + cosAlpha * poly_->points[0].x() - sinAlpha * poly_->points[0].y() + direction_[0], + sinAlpha * poly_->points[0].x() + cosAlpha * poly_->points[0].y() + direction_[1], + direction_[2] + } + }); + + // features connecting the top and bottom + if (alpha_ == 0) { + for (const auto & pt: poly_->points) { + std::vector> line = { + {pt.x(), pt.y(), 0.0}, + {pt.x() + direction_[0], pt.y() + direction_[1], direction_[2]} + }; + features.push_back(line); + } + } else { + // Alright, we need to chop the lines on which the polygon corners are + // sitting into pieces. How long? About edge_size. For the starting point + // (x0, y0, z0) height h and angle alpha, the lines are given by + // + // f(beta) = ( + // cos(alpha*beta) x0 - sin(alpha*beta) y0, + // sin(alpha*beta) x0 + cos(alpha*beta) y0, + // z0 + beta * h + // ) + // + // with beta in [0, 1]. The length from beta0 till beta1 is then + // + // l = sqrt(alpha^2 (x0^2 + y0^2) + h^2) * (beta1 - beta0). + // + const double height = direction_[2]; + for (const auto & pt: poly_->points) { + const double l = sqrt(alpha_*alpha_ * (pt.x()*pt.x() + pt.y()*pt.y()) + height*height); + assert(edge_size_ > 0.0); + const size_t n = int(l / edge_size_ - 0.5) + 1; + std::vector> line = { + {pt.x(), pt.y(), 0.0}, + }; + for (size_t i=0; i < n; i++) { + const double beta = double(i+1) / n; + const double sinAB = sin(alpha_*beta); + const double cosAB = cos(alpha_*beta); + line.push_back({ + cosAB * pt.x() - sinAB * pt.y(), + sinAB * pt.x() + cosAB * pt.y(), + beta * height + }); + } + features.push_back(line); + } + } + + return features; + }; + + private: + const std::shared_ptr poly_; + const std::array direction_; + const double alpha_; + const double edge_size_; +}; + + +class ring_extrude: public pygalmesh::DomainBase { + public: + ring_extrude( + const std::shared_ptr & poly, + const double edge_size + ): + poly_(poly), + edge_size_(edge_size) + { + assert(edge_size > 0.0); + } + + virtual ~ring_extrude() = default; + + virtual + double + eval(const std::array & x) const + { + const double r = sqrt(x[0]*x[0] + x[1]*x[1]); + const double z = x[2]; + + return poly_->is_inside({r, z}) ? -1.0 : 1.0; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + double max = 0.0; + for (const auto & pt: poly_->points) { + const double nrm1 = pt.x()*pt.x() + pt.y()*pt.y(); + if (nrm1 > max) { + max = nrm1; + } + } + return max; + } + + virtual + std::vector>> + get_features() const + { + std::vector>> features = {}; + + for (const auto & pt: poly_->points) { + const double r = pt.x(); + const double circ = 2 * 3.14159265359 * r; + const size_t n = int(circ / edge_size_ - 0.5) + 1; + std::vector> line; + for (size_t i=0; i < n; i++) { + const double alpha = (2 * 3.14159265359 * i) / n; + line.push_back({ + r * cos(alpha), + r * sin(alpha), + pt.y() + }); + } + line.push_back(line.front()); + features.push_back(line); + } + return features; + } + + private: + const std::shared_ptr poly_; + const double edge_size_; +}; + +} // namespace pygalmesh + +#endif // POLYGON2D_HPP diff --git a/src/primitives.hpp b/src/primitives.hpp new file mode 100644 index 0000000..f3e66c3 --- /dev/null +++ b/src/primitives.hpp @@ -0,0 +1,453 @@ +#ifndef PRIMITIVES_HPP +#define PRIMITIVES_HPP + +#include "domain.hpp" + +#include +#include + +namespace pygalmesh { + +class Ball: public pygalmesh::DomainBase +{ + public: + Ball( + const std::array & x0, + const double radius + ): + x0_(x0), + radius_(radius) + { + assert(x0_.size() == 3); + } + + virtual ~Ball() = default; + + virtual + double + eval(const std::array & x) const + { + const double xx0 = x[0] - x0_[0]; + const double yy0 = x[1] - x0_[1]; + const double zz0 = x[2] - x0_[2]; + return xx0*xx0 + yy0*yy0 + zz0*zz0 - radius_*radius_; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + const double x0_nrm = sqrt(x0_[0]*x0_[0] + x0_[1]*x0_[1] + x0_[2]*x0_[2]); + return (x0_nrm + radius_) * (x0_nrm + radius_); + } + + private: + const std::array x0_; + const double radius_; +}; + + +class Cuboid: public pygalmesh::DomainBase +{ + public: + Cuboid( + const std::array & x0, + const std::array & x1 + ): + x0_(x0), + x1_(x1) + { + } + + virtual ~Cuboid() = default; + + virtual + double + eval(const std::array & x) const + { + // TODO differentiable expression? + return std::max(std::max( + (x[0] - x0_[0]) * (x[0] - x1_[0]), + (x[1] - x0_[1]) * (x[1] - x1_[1]) + ), + (x[2] - x0_[2]) * (x[2] - x1_[2]) + ); + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + const double x0_nrm2 = x0_[0]*x0_[0] + x0_[1]*x0_[1] + x0_[2]*x0_[2]; + const double x1_nrm2 = x1_[0]*x1_[0] + x1_[1]*x1_[1] + x1_[2]*x1_[2]; + return std::max({x0_nrm2, x1_nrm2}); + } + + virtual + std::vector>> + get_features() const + { + std::vector> corners = { + {x0_[0], x0_[1], x0_[2]}, + {x1_[0], x0_[1], x0_[2]}, + {x0_[0], x1_[1], x0_[2]}, + {x0_[0], x0_[1], x1_[2]}, + {x1_[0], x1_[1], x0_[2]}, + {x1_[0], x0_[1], x1_[2]}, + {x0_[0], x1_[1], x1_[2]}, + {x1_[0], x1_[1], x1_[2]} + }; + return { + {corners[0], corners[1]}, + {corners[0], corners[2]}, + {corners[0], corners[3]}, + {corners[1], corners[4]}, + {corners[1], corners[5]}, + {corners[2], corners[4]}, + {corners[2], corners[6]}, + {corners[3], corners[5]}, + {corners[3], corners[6]}, + {corners[4], corners[7]}, + {corners[5], corners[7]}, + {corners[6], corners[7]} + }; + }; + + private: + const std::array x0_; + const std::array x1_; +}; + + +class Ellipsoid: public pygalmesh::DomainBase +{ + public: + Ellipsoid( + const std::array & x0, + const double a0, + const double a1, + const double a2 + ): + x0_(x0), + a0_2_(a0*a0), + a1_2_(a1*a1), + a2_2_(a2*a1) + { + } + + virtual ~Ellipsoid() = default; + + virtual + double + eval(const std::array & x) const + { + const double xx0 = x[0] - x0_[0]; + const double yy0 = x[1] - x0_[1]; + const double zz0 = x[2] - x0_[2]; + return xx0*xx0/a0_2_ + yy0*yy0/a1_2_ + zz0*zz0/a2_2_ - 1.0; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + return std::max({a0_2_, a1_2_, a2_2_}); + } + + private: + const std::array x0_; + const double a0_2_; + const double a1_2_; + const double a2_2_; +}; + + +class Cylinder: public pygalmesh::DomainBase +{ + public: + Cylinder( + const double z0, + const double z1, + const double radius, + const double feature_edge_length + ): + z0_(z0), + z1_(z1), + radius_(radius), + feature_edge_length_(feature_edge_length) + { + assert(z1_ > z0_); + } + + virtual ~Cylinder() = default; + + virtual + double + eval(const std::array & x) const + { + return (z0_ < x[2] && x[2] < z1_) ? + x[0]*x[0] + x[1]*x[1] - radius_*radius_ : + 1.0; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + const double zmax = std::max({abs(z0_), abs(z1_)}); + return zmax*zmax + radius_*radius_; + } + + virtual + std::vector>> + get_features() const + { + const double pi = 3.1415926535897932384; + const size_t n = 2 * pi * radius_ / feature_edge_length_; + std::vector> circ0(n+1); + std::vector> circ1(n+1); + for (size_t i=0; i < n; i++) { + const double c = radius_ * cos((2*pi * i) / n); + const double s = radius_ * sin((2*pi * i) / n); + circ0[i] = {c, s, z0_}; + circ1[i] = {c, s, z1_}; + } + // close the circles + circ0[n] = circ0[0]; + circ1[n] = circ1[0]; + return {circ0, circ1}; + }; + + private: + const double z0_; + const double z1_; + const double radius_; + const double feature_edge_length_; +}; + + +class Cone: public pygalmesh::DomainBase +{ + public: + Cone( + const double radius, + const double height, + const double feature_edge_length + ): + radius_(radius), + height_(height), + feature_edge_length_(feature_edge_length) + { + assert(radius_ > 0.0); + assert(height_ > 0.0); + } + + virtual ~Cone() = default; + + virtual + double + eval(const std::array & x) const + { + const double rad = radius_ * (1.0 - x[2] / height_); + + return (0.0 < x[2] && x[2] < height_) ? + x[0]*x[0] + x[1]*x[1] - rad*rad : + 1.0; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + const double max = std::max({radius_, height_}); + return max*max; + } + + virtual + std::vector>> + get_features() const + { + const double pi = 3.1415926535897932384; + const size_t n = 2 * pi * radius_ / feature_edge_length_; + std::vector> circ0(n+1); + for (size_t i=0; i < n; i++) { + const double c = radius_ * cos((2*pi * i) / n); + const double s = radius_ * sin((2*pi * i) / n); + circ0[i] = {c, s, 0.0}; + } + circ0[n] = circ0[0]; + return {circ0}; + }; + + private: + const double radius_; + const double height_; + const double feature_edge_length_; +}; + + +class Tetrahedron: public pygalmesh::DomainBase +{ + public: + Tetrahedron( + const std::array & x0, + const std::array & x1, + const std::array & x2, + const std::array & x3 + ): + x0_(Eigen::Vector3d(x0[0], x0[1], x0[2])), + x1_(Eigen::Vector3d(x1[0], x1[1], x1[2])), + x2_(Eigen::Vector3d(x2[0], x2[1], x2[2])), + x3_(Eigen::Vector3d(x3[0], x3[1], x3[2])) + { + } + + virtual ~Tetrahedron() = default; + + bool isOnSameSide( + const Eigen::Vector3d & v0, + const Eigen::Vector3d & v1, + const Eigen::Vector3d & v2, + const Eigen::Vector3d & v3, + const Eigen::Vector3d & p + ) const + { + const auto normal = (v1 - v0).cross(v2 - v0); + const double dot_v3 = normal.dot(v3 - v0); + const double dot_p = normal.dot(p - v0); + return ( + (dot_v3 > 0 && dot_p > 0) || (dot_v3 < 0 && dot_p < 0) + ); + } + + virtual + double + eval(const std::array & x) const + { + // TODO continuous expression + Eigen::Vector3d pvec(x.data()); + const bool a = + isOnSameSide(x0_, x1_, x2_, x3_, pvec) && + isOnSameSide(x1_, x2_, x3_, x0_, pvec) && + isOnSameSide(x2_, x3_, x0_, x1_, pvec) && + isOnSameSide(x3_, x0_, x1_, x2_, pvec); + return a ? -1.0 : 1.0; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + return std::max({ + x0_.dot(x0_), + x1_.dot(x1_), + x2_.dot(x2_), + x3_.dot(x3_) + }); + } + + virtual + std::vector>> + get_features() const + { + std::vector> pts = { + {x0_[0], x0_[1], x0_[2]}, + {x1_[0], x1_[1], x1_[2]}, + {x2_[0], x2_[1], x2_[2]}, + {x3_[0], x3_[1], x3_[2]} + }; + return { + {pts[0], pts[1]}, + {pts[0], pts[2]}, + {pts[0], pts[3]}, + {pts[1], pts[2]}, + {pts[1], pts[3]}, + {pts[2], pts[3]} + }; + }; + + private: + const Eigen::Vector3d x0_; + const Eigen::Vector3d x1_; + const Eigen::Vector3d x2_; + const Eigen::Vector3d x3_; +}; + + +class Torus: public pygalmesh::DomainBase +{ + public: + Torus( + const double major_radius, + const double minor_radius + ): + major_radius_(major_radius), + minor_radius_(minor_radius) + { + } + + virtual ~Torus() = default; + + virtual + double + eval(const std::array & x) const + { + const double r = sqrt(x[0]*x[0] + x[1]*x[1]); + return ( + (r - major_radius_)*(r - major_radius_) + x[2]*x[2] + - minor_radius_*minor_radius_ + ); + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + return (major_radius_ + minor_radius_)*(major_radius_ + minor_radius_); + } + + private: + const double major_radius_; + const double minor_radius_; +}; + + +class HalfSpace: public pygalmesh::DomainBase +{ + public: + HalfSpace( + const std::array & n, + const double alpha, + const double bounding_sphere_squared_radius + ): + n_(n), + alpha_(alpha), + bounding_sphere_squared_radius_(bounding_sphere_squared_radius) + { + } + + virtual ~HalfSpace() = default; + + virtual + double + eval(const std::array & x) const + { + return n_[0]*x[0] + n_[1]*x[1] + n_[2]*x[2] - alpha_; + } + + virtual + double + get_bounding_sphere_squared_radius() const + { + return bounding_sphere_squared_radius_; + } + + private: + const std::array n_; + const double alpha_; + const double bounding_sphere_squared_radius_; +}; + +} // namespace pygalmesh + +#endif // PRIMITIVES_HPP diff --git a/src/pybind11.cpp b/src/pybind11.cpp new file mode 100644 index 0000000..0ade4be --- /dev/null +++ b/src/pybind11.cpp @@ -0,0 +1,265 @@ +#include "domain.hpp" +#include "generate.hpp" +#include "generate_from_off.hpp" +#include "generate_surface_mesh.hpp" +#include "polygon2d.hpp" +#include "primitives.hpp" + +#include +#include + +namespace py = pybind11; +using namespace pygalmesh; + + +// https://pybind11.readthedocs.io/en/stable/advanced/classes.html#overriding-virtual-functions-in-python +class PyDomainBase: public DomainBase { +public: + using DomainBase::DomainBase; + + double + eval(const std::array & x) const override { + PYBIND11_OVERLOAD_PURE(double, DomainBase, eval, x); + } + + double + get_bounding_sphere_squared_radius() const override { + PYBIND11_OVERLOAD_PURE(double, DomainBase, get_bounding_sphere_squared_radius); + } + + // std::vector>> + // get_features() const override { + // PYBIND11_OVERLOAD( + // std::vector>>, + // DomainBase, + // get_features, + // 0.0 + // ); + // } +}; + + +PYBIND11_MODULE(_pygalmesh, m) { + // m.doc() = "documentation string"; + + // Domain base. + // shared_ptr b/c of + // + py::class_>(m, "DomainBase") + .def(py::init<>()) + .def("eval", &DomainBase::eval) + .def("get_bounding_sphere_squared_radius", &DomainBase::get_bounding_sphere_squared_radius) + .def("get_features", &DomainBase::get_features); + + // Domain transformations + py::class_>(m, "Translate") + .def(py::init< + const std::shared_ptr &, + const std::array & + >()) + .def("eval", &Translate::eval) + .def("translate_features", &Translate::translate_features) + .def("get_bounding_sphere_squared_radius", &Translate::get_bounding_sphere_squared_radius) + .def("get_features", &Translate::get_features); + + py::class_>(m, "Rotate") + .def(py::init< + const std::shared_ptr &, + const std::array &, + const double + >()) + .def("eval", &Rotate::eval) + .def("rotate", &Rotate::rotate) + .def("rotate_features", &Rotate::rotate_features) + .def("get_bounding_sphere_squared_radius", &Rotate::get_bounding_sphere_squared_radius) + .def("get_features", &Rotate::get_features); + + py::class_>(m, "Scale") + .def(py::init< + std::shared_ptr &, + const double + >()) + .def("eval", &Scale::eval) + .def("scale_features", &Scale::scale_features) + .def("get_bounding_sphere_squared_radius", &Scale::get_bounding_sphere_squared_radius) + .def("get_features", &Scale::get_features); + + py::class_>(m, "Stretch") + .def(py::init< + std::shared_ptr &, + const std::array & + >()) + .def("eval", &Stretch::eval) + .def("stretch_features", &Stretch::stretch_features) + .def("get_bounding_sphere_squared_radius", &Stretch::get_bounding_sphere_squared_radius) + .def("get_features", &Stretch::get_features); + + py::class_>(m, "Intersection") + .def(py::init< + std::vector> & + >()) + .def("eval", &Intersection::eval) + .def("get_bounding_sphere_squared_radius", &Intersection::get_bounding_sphere_squared_radius) + .def("get_features", &Intersection::get_features); + + py::class_>(m, "Union") + .def(py::init< + std::vector> & + >()) + .def("eval", &Union::eval) + .def("get_bounding_sphere_squared_radius", &Union::get_bounding_sphere_squared_radius) + .def("get_features", &Union::get_features); + + py::class_>(m, "Difference") + .def(py::init< + std::shared_ptr &, + std::shared_ptr & + >()) + .def("eval", &Difference::eval) + .def("get_bounding_sphere_squared_radius", &Difference::get_bounding_sphere_squared_radius) + .def("get_features", &Difference::get_features); + + // Primitives + py::class_>(m, "Ball") + .def(py::init< + const std::array &, + const double + >()) + .def("eval", &Ball::eval) + .def("get_bounding_sphere_squared_radius", &Ball::get_bounding_sphere_squared_radius); + + py::class_>(m, "Cuboid") + .def(py::init< + const std::array &, + const std::array & + >()) + .def("eval", &Cuboid::eval) + .def("get_bounding_sphere_squared_radius", &Cuboid::get_bounding_sphere_squared_radius) + .def("get_features", &Cuboid::get_features); + + py::class_>(m, "Ellipsoid") + .def(py::init< + const std::array &, + const double, + const double, + const double + >()) + .def("eval", &Ellipsoid::eval) + .def("get_bounding_sphere_squared_radius", &Ellipsoid::get_bounding_sphere_squared_radius) + .def("get_features", &Ellipsoid::get_features); + + py::class_>(m, "Cylinder") + .def(py::init< + const double, + const double, + const double, + const double + >()) + .def("eval", &Cylinder::eval) + .def("get_bounding_sphere_squared_radius", &Cylinder::get_bounding_sphere_squared_radius) + .def("get_features", &Cylinder::get_features); + + py::class_>(m, "Cone") + .def(py::init< + const double, + const double, + const double + >()) + .def("eval", &Cone::eval) + .def("get_bounding_sphere_squared_radius", &Cone::get_bounding_sphere_squared_radius) + .def("get_features", &Cone::get_features); + + py::class_>(m, "Tetrahedron") + .def(py::init< + const std::array &, + const std::array &, + const std::array &, + const std::array & + >()) + .def("eval", &Tetrahedron::eval) + .def("get_bounding_sphere_squared_radius", &Tetrahedron::get_bounding_sphere_squared_radius) + .def("get_features", &Tetrahedron::get_features); + + py::class_>(m, "Torus") + .def(py::init< + const double, + const double + >()) + .def("eval", &Torus::eval) + .def("get_bounding_sphere_squared_radius", &Torus::get_bounding_sphere_squared_radius) + .def("get_features", &Torus::get_features); + + py::class_>(m, "HalfSpace") + .def(py::init< + const std::array &, + const double, + const double + >()) + .def("eval", &HalfSpace::eval) + .def("get_bounding_sphere_squared_radius", &HalfSpace::get_bounding_sphere_squared_radius); + + // polygon2d + py::class_>(m, "Polygon2D") + .def(py::init< + const std::vector> & + >()) + .def("vector_to_cgal_points", &Polygon2D::vector_to_cgal_points) + .def("is_inside", &Polygon2D::is_inside); + + py::class_>(m, "Extrude") + .def(py::init< + const std::shared_ptr &, + const std::array &, + const double, + const double + >(), + py::arg("poly"), + py::arg("direction"), + py::arg("alpha") = 0.0, + py::arg("edge_size") = 0.0 + ) + .def("eval", &Extrude::eval) + .def("get_bounding_sphere_squared_radius", &Extrude::get_bounding_sphere_squared_radius) + .def("get_features", &Extrude::get_features); + + py::class_>(m, "RingExtrude") + .def(py::init< + const std::shared_ptr &, + const double + >()) + .def("eval", &ring_extrude::eval) + .def("get_bounding_sphere_squared_radius", &ring_extrude::get_bounding_sphere_squared_radius) + .def("get_features", &ring_extrude::get_features); + + // functions + m.def("generate_from_off", &generate_from_off); + m.def( + "generate_mesh", &generate_mesh, + py::arg("domain"), + py::arg("outfile"), + py::arg("feature_edges") = std::vector>>(), + py::arg("bounding_sphere_radius") = 0.0, + py::arg("boundary_precision") = 1.0e-4, + py::arg("lloyd") = false, + py::arg("odt") = false, + py::arg("perturb") = true, + py::arg("exude") = true, + py::arg("edge_size") = 0.0, + py::arg("facet_angle") = 0.0, + py::arg("facet_size") = 0.0, + py::arg("facet_distance") = 0.0, + py::arg("cell_radius_edge_ratio") = 0.0, + py::arg("cell_size") = 0.0, + py::arg("verbose") = true + ); + m.def( + "generate_surface_mesh", &generate_surface_mesh, + py::arg("domain"), + py::arg("outfile"), + py::arg("bounding_sphere_radius") = 0.0, + py::arg("angle_bound") = 0.0, + py::arg("radius_bound") = 0.0, + py::arg("distance_bound") = 0.0, + py::arg("verbose") = true + ); +} diff --git a/test/test_pygalmesh.py b/test/test_pygalmesh.py new file mode 100644 index 0000000..977840a --- /dev/null +++ b/test/test_pygalmesh.py @@ -0,0 +1,680 @@ +# -*- coding: utf-8 -*- +# +import numpy +import meshio + +# pylint: disable=import-error +import pygalmesh + + +def _row_dot(a, b): + # http://stackoverflow.com/a/26168677/353337 + return numpy.einsum("ij, ij->i", a, b) + + +def compute_volumes(vertices, tets): + cell_coords = vertices[tets] + + a = cell_coords[:, 1, :] - cell_coords[:, 0, :] + b = cell_coords[:, 2, :] - cell_coords[:, 0, :] + c = cell_coords[:, 3, :] - cell_coords[:, 0, :] + + # omega = + omega = _row_dot(a, numpy.cross(b, c)) + + # https://en.wikipedia.org/wiki/Tetrahedron#Volume + return abs(omega) / 6.0 + + +def compute_triangle_areas(vertices, triangles): + e0 = vertices[triangles[:, 1]] - vertices[triangles[:, 0]] + e1 = vertices[triangles[:, 2]] - vertices[triangles[:, 1]] + + e0_cross_e1 = numpy.cross(e0, e1) + areas = 0.5 * numpy.sqrt(_row_dot(e0_cross_e1, e0_cross_e1)) + + return areas + + +def test_ball(): + s = pygalmesh.Ball([0.0, 0.0, 0.0], 1.0) + pygalmesh.generate_mesh(s, "out.mesh", cell_size=0.2, verbose=False) + + mesh = meshio.read("out.mesh") + + assert abs(max(mesh.points[:, 0]) - 1.0) < 0.02 + assert abs(min(mesh.points[:, 0]) + 1.0) < 0.02 + assert abs(max(mesh.points[:, 1]) - 1.0) < 0.02 + assert abs(min(mesh.points[:, 1]) + 1.0) < 0.02 + assert abs(max(mesh.points[:, 2]) - 1.0) < 0.02 + assert abs(min(mesh.points[:, 2]) + 1.0) < 0.02 + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 4.0 / 3.0 * numpy.pi) < 0.15 + return + + +def test_balls_union(): + radius = 1.0 + displacement = 0.5 + s0 = pygalmesh.Ball([displacement, 0, 0], radius) + s1 = pygalmesh.Ball([-displacement, 0, 0], radius) + u = pygalmesh.Union([s0, s1]) + + a = numpy.sqrt(radius ** 2 - displacement ** 2) + edge_size = 0.1 + n = int(2 * numpy.pi * a / edge_size) + circ = [ + [0.0, a * numpy.cos(i * 2 * numpy.pi / n), a * numpy.sin(i * 2 * numpy.pi / n)] + for i in range(n) + ] + circ.append(circ[0]) + + pygalmesh.generate_mesh( + u, + "out.mesh", + feature_edges=[circ], + cell_size=0.15, + edge_size=edge_size, + verbose=False, + ) + + mesh = meshio.read("out.mesh") + + assert abs(max(mesh.points[:, 0]) - (radius + displacement)) < 0.02 + assert abs(min(mesh.points[:, 0]) + (radius + displacement)) < 0.02 + assert abs(max(mesh.points[:, 1]) - radius) < 0.02 + assert abs(min(mesh.points[:, 1]) + radius) < 0.02 + assert abs(max(mesh.points[:, 2]) - radius) < 0.02 + assert abs(min(mesh.points[:, 2]) + radius) < 0.02 + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + h = radius - displacement + ref_vol = 2 * ( + 4.0 / 3.0 * numpy.pi * radius ** 3 - h * numpy.pi / 6.0 * (3 * a ** 2 + h ** 2) + ) + + assert abs(vol - ref_vol) < 0.1 + + return + + +def test_balls_intersection(): + radius = 1.0 + displacement = 0.5 + s0 = pygalmesh.Ball([displacement, 0, 0], radius) + s1 = pygalmesh.Ball([-displacement, 0, 0], radius) + u = pygalmesh.Intersection([s0, s1]) + + a = numpy.sqrt(radius ** 2 - displacement ** 2) + edge_size = 0.1 + n = int(2 * numpy.pi * a / edge_size) + circ = [ + [0.0, a * numpy.cos(i * 2 * numpy.pi / n), a * numpy.sin(i * 2 * numpy.pi / n)] + for i in range(n) + ] + circ.append(circ[0]) + + pygalmesh.generate_mesh( + u, + "out.mesh", + feature_edges=[circ], + cell_size=0.15, + edge_size=edge_size, + verbose=False, + ) + + mesh = meshio.read("out.mesh") + + assert abs(max(mesh.points[:, 0]) - (radius - displacement)) < 0.02 + assert abs(min(mesh.points[:, 0]) + (radius - displacement)) < 0.02 + assert abs(max(mesh.points[:, 1]) - a) < 0.02 + assert abs(min(mesh.points[:, 1]) + a) < 0.02 + assert abs(max(mesh.points[:, 2]) - a) < 0.02 + assert abs(min(mesh.points[:, 2]) + a) < 0.02 + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + h = radius - displacement + ref_vol = 2 * (h * numpy.pi / 6.0 * (3 * a ** 2 + h ** 2)) + + assert abs(vol - ref_vol) < 0.1 + + return + + +# pylint: disable=too-many-locals +def test_balls_difference(): + radius = 1.0 + displacement = 0.5 + s0 = pygalmesh.Ball([displacement, 0, 0], radius) + s1 = pygalmesh.Ball([-displacement, 0, 0], radius) + u = pygalmesh.Difference(s0, s1) + + a = numpy.sqrt(radius ** 2 - displacement ** 2) + edge_size = 0.15 + n = int(2 * numpy.pi * a / edge_size) + circ = [ + [0.0, a * numpy.cos(i * 2 * numpy.pi / n), a * numpy.sin(i * 2 * numpy.pi / n)] + for i in range(n) + ] + circ.append(circ[0]) + + pygalmesh.generate_mesh( + u, + "out.mesh", + feature_edges=[circ], + cell_size=0.15, + edge_size=edge_size, + facet_angle=25, + facet_size=0.15, + cell_radius_edge_ratio=2.0, + verbose=False, + ) + + mesh = meshio.read("out.mesh") + + tol = 0.02 + assert abs(max(mesh.points[:, 0]) - (radius + displacement)) < tol + assert abs(min(mesh.points[:, 0]) - 0.0) < tol + assert abs(max(mesh.points[:, 1]) - radius) < tol + assert abs(min(mesh.points[:, 1]) + radius) < tol + assert abs(max(mesh.points[:, 2]) - radius) < tol + assert abs(min(mesh.points[:, 2]) + radius) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + h = radius - displacement + ref_vol = 4.0 / 3.0 * numpy.pi * radius ** 3 - 2 * h * numpy.pi / 6.0 * ( + 3 * a ** 2 + h ** 2 + ) + + assert abs(vol - ref_vol) < 0.05 + + return + + +def test_cuboids_intersection(): + c0 = pygalmesh.Cuboid([0, 0, -0.5], [3, 3, 0.5]) + c1 = pygalmesh.Cuboid([1, 1, -2], [2, 2, 2]) + u = pygalmesh.Intersection([c0, c1]) + + # In CGAL, feature edges must not intersect, and that's a problem here: The + # intersection edges of the cuboids share eight points with the edges of + # the tall and skinny cuboid. + # eps = 1.0e-2 + # extra_features = [ + # [[1.0, 1.0 + eps, 0.5], [1.0, 2.0 - eps, 0.5]], + # [[1.0 + eps, 2.0, 0.5], [2.0 - eps, 2.0, 0.5]], + # [[2.0, 2.0 - eps, 0.5], [2.0, 1.0 + eps, 0.5]], + # [[2.0 - eps, 1.0, 0.5], [1.0 + eps, 1.0, 0.5]], + # ] + + pygalmesh.generate_mesh(u, "out.mesh", cell_size=0.1, edge_size=0.1, verbose=False) + + mesh = meshio.read("out.mesh") + + # filter the vertices that belong to cells + verts = mesh.points[numpy.unique(mesh.cells["tetra"])] + + tol = 1.0e-2 + assert abs(max(verts[:, 0]) - 2.0) < tol + assert abs(min(verts[:, 0]) - 1.0) < tol + assert abs(max(verts[:, 1]) - 2.0) < tol + assert abs(min(verts[:, 1]) - 1.0) < tol + assert abs(max(verts[:, 2]) - 0.5) < 0.05 + assert abs(min(verts[:, 2]) + 0.5) < 0.05 + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 1.0) < 0.05 + + return + + +def test_cuboids_union(): + c0 = pygalmesh.Cuboid([0, 0, -0.5], [3, 3, 0.5]) + c1 = pygalmesh.Cuboid([1, 1, -2], [2, 2, 2]) + u = pygalmesh.Union([c0, c1]) + + pygalmesh.generate_mesh(u, "out.mesh", cell_size=0.2, edge_size=0.2, verbose=False) + + mesh = meshio.read("out.mesh") + + # filter the vertices that belong to cells + verts = mesh.points[numpy.unique(mesh.cells["tetra"])] + + tol = 1.0e-2 + assert abs(max(verts[:, 0]) - 3.0) < tol + assert abs(min(verts[:, 0]) - 0.0) < tol + assert abs(max(verts[:, 1]) - 3.0) < tol + assert abs(min(verts[:, 1]) - 0.0) < tol + assert abs(max(verts[:, 2]) - 2.0) < tol + assert abs(min(verts[:, 2]) + 2.0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 12.0) < 0.1 + return + + +def test_cuboid(): + s0 = pygalmesh.Cuboid([0, 0, 0], [1, 2, 3]) + pygalmesh.generate_mesh(s0, "out.mesh", edge_size=0.1, verbose=False) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-3 + assert abs(max(mesh.points[:, 0]) - 1.0) < tol + assert abs(min(mesh.points[:, 0]) + 0.0) < tol + assert abs(max(mesh.points[:, 1]) - 2.0) < tol + assert abs(min(mesh.points[:, 1]) + 0.0) < tol + assert abs(max(mesh.points[:, 2]) - 3.0) < tol + assert abs(min(mesh.points[:, 2]) + 0.0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 6.0) < tol + return + + +def test_cone(): + base_radius = 1.0 + height = 2.0 + edge_size = 0.1 + s0 = pygalmesh.Cone(base_radius, height, edge_size) + pygalmesh.generate_mesh( + s0, "out.mesh", cell_size=0.1, edge_size=edge_size, verbose=False + ) + + mesh = meshio.read("out.mesh") + + tol = 2.0e-1 + assert abs(max(mesh.points[:, 0]) - base_radius) < tol + assert abs(min(mesh.points[:, 0]) + base_radius) < tol + assert abs(max(mesh.points[:, 1]) - base_radius) < tol + assert abs(min(mesh.points[:, 1]) + base_radius) < tol + assert abs(max(mesh.points[:, 2]) - height) < tol + assert abs(min(mesh.points[:, 2]) + 0.0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + ref_vol = numpy.pi * base_radius * base_radius / 3.0 * height + assert abs(vol - ref_vol) < tol + return + + +def test_cylinder(): + radius = 1.0 + z0 = 0.0 + z1 = 1.0 + edge_length = 0.1 + s0 = pygalmesh.Cylinder(z0, z1, radius, edge_length) + pygalmesh.generate_mesh( + s0, "out.mesh", cell_size=0.1, edge_size=edge_length, verbose=False + ) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-1 + assert abs(max(mesh.points[:, 0]) - radius) < tol + assert abs(min(mesh.points[:, 0]) + radius) < tol + assert abs(max(mesh.points[:, 1]) - radius) < tol + assert abs(min(mesh.points[:, 1]) + radius) < tol + assert abs(max(mesh.points[:, 2]) - z1) < tol + assert abs(min(mesh.points[:, 2]) + z0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + ref_vol = numpy.pi * radius * radius * (z1 - z0) + assert abs(vol - ref_vol) < tol + return + + +def test_tetrahedron(): + s0 = pygalmesh.Tetrahedron( + [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0] + ) + pygalmesh.generate_mesh(s0, "out.mesh", cell_size=0.1, edge_size=0.1, verbose=False) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-4 + assert abs(max(mesh.points[:, 0]) - 1.0) < tol + assert abs(min(mesh.points[:, 0]) + 0.0) < tol + assert abs(max(mesh.points[:, 1]) - 1.0) < tol + assert abs(min(mesh.points[:, 1]) + 0.0) < tol + assert abs(max(mesh.points[:, 2]) - 1.0) < tol + assert abs(min(mesh.points[:, 2]) + 0.0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 1.0 / 6.0) < tol + return + + +def test_torus(): + major_radius = 1.0 + minor_radius = 0.5 + s0 = pygalmesh.Torus(major_radius, minor_radius) + pygalmesh.generate_mesh(s0, "out.mesh", cell_size=0.1, verbose=False) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-2 + radii_sum = major_radius + minor_radius + assert abs(max(mesh.points[:, 0]) - radii_sum) < tol + assert abs(min(mesh.points[:, 0]) + radii_sum) < tol + assert abs(max(mesh.points[:, 1]) - radii_sum) < tol + assert abs(min(mesh.points[:, 1]) + radii_sum) < tol + assert abs(max(mesh.points[:, 2]) - minor_radius) < tol + assert abs(min(mesh.points[:, 2]) + minor_radius) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + ref_vol = (numpy.pi * minor_radius * minor_radius) * (2 * numpy.pi * major_radius) + assert abs(vol - ref_vol) < 1.0e-1 + return + + +def test_custom_function(): + class Hyperboloid(pygalmesh.DomainBase): + def __init__(self, edge_size): + super(Hyperboloid, self).__init__() + self.z0 = -1.0 + self.z1 = 1.0 + self.waist_radius = 0.5 + self.edge_size = edge_size + return + + def eval(self, x): + if self.z0 < x[2] and x[2] < self.z1: + return x[0] ** 2 + x[1] ** 2 - (x[2] ** 2 + self.waist_radius) ** 2 + return 1.0 + + def get_bounding_sphere_squared_radius(self): + z_max = max(abs(self.z0), abs(self.z1)) + r_max = z_max ** 2 + self.waist_radius + return r_max * r_max + z_max * z_max + + def get_features(self): + radius0 = self.z0 ** 2 + self.waist_radius + n0 = int(2 * numpy.pi * radius0 / self.edge_size) + circ0 = [ + [ + radius0 * numpy.cos((2 * numpy.pi * k) / n0), + radius0 * numpy.sin((2 * numpy.pi * k) / n0), + self.z0, + ] + for k in range(n0) + ] + circ0.append(circ0[0]) + + radius1 = self.z1 ** 2 + self.waist_radius + n1 = int(2 * numpy.pi * radius1 / self.edge_size) + circ1 = [ + [ + radius1 * numpy.cos((2 * numpy.pi * k) / n1), + radius1 * numpy.sin((2 * numpy.pi * k) / n1), + self.z1, + ] + for k in range(n1) + ] + circ1.append(circ1[0]) + return [circ0, circ1] + + edge_size = 0.12 + d = Hyperboloid(edge_size) + + pygalmesh.generate_mesh( + d, "out.mesh", cell_size=0.1, edge_size=edge_size, verbose=False + ) + + mesh = meshio.read("out.mesh") + + # TODO check the reference values + tol = 1.0e-1 + assert abs(max(mesh.points[:, 0]) - 1.4) < tol + assert abs(min(mesh.points[:, 0]) + 1.4) < tol + assert abs(max(mesh.points[:, 1]) - 1.4) < tol + assert abs(min(mesh.points[:, 1]) + 1.4) < tol + assert abs(max(mesh.points[:, 2]) - 1.0) < tol + assert abs(min(mesh.points[:, 2]) + 1.0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 2 * numpy.pi * 47.0 / 60.0) < 0.15 + return + + +def test_scaling(): + alpha = 1.3 + s = pygalmesh.Scale(pygalmesh.Cuboid([0, 0, 0], [1, 2, 3]), alpha) + pygalmesh.generate_mesh(s, "out.mesh", cell_size=0.2, edge_size=0.1, verbose=False) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-3 + assert abs(max(mesh.points[:, 0]) - 1 * alpha) < tol + assert abs(min(mesh.points[:, 0]) + 0.0) < tol + assert abs(max(mesh.points[:, 1]) - 2 * alpha) < tol + assert abs(min(mesh.points[:, 1]) + 0.0) < tol + assert abs(max(mesh.points[:, 2]) - 3 * alpha) < tol + assert abs(min(mesh.points[:, 2]) + 0.0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 6.0 * alpha ** 3) < tol + return + + +def test_stretch(): + alpha = 2.0 + s = pygalmesh.Stretch(pygalmesh.Cuboid([0, 0, 0], [1, 2, 3]), [alpha, 0.0, 0.0]) + pygalmesh.generate_mesh(s, "out.mesh", cell_size=0.2, edge_size=0.2, verbose=False) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-3 + assert abs(max(mesh.points[:, 0]) - alpha) < tol + assert abs(min(mesh.points[:, 0]) + 0.0) < tol + assert abs(max(mesh.points[:, 1]) - 2.0) < tol + assert abs(min(mesh.points[:, 1]) + 0.0) < tol + assert abs(max(mesh.points[:, 2]) - 3.0) < tol + assert abs(min(mesh.points[:, 2]) + 0.0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 12.0) < tol + + return + + +def test_rotation(): + s0 = pygalmesh.Rotate( + pygalmesh.Cuboid([0, 0, 0], [1, 2, 3]), [1.0, 0.0, 0.0], numpy.pi / 12.0 + ) + pygalmesh.generate_mesh(s0, "out.mesh", cell_size=0.1, edge_size=0.1, verbose=False) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-3 + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 6.0) < tol + return + + +def test_translation(): + s0 = pygalmesh.Translate(pygalmesh.Cuboid([0, 0, 0], [1, 2, 3]), [1.0, 0.0, 0.0]) + pygalmesh.generate_mesh(s0, "out.mesh", cell_size=0.1, edge_size=0.1, verbose=False) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-3 + assert abs(max(mesh.points[:, 0]) - 2.0) < tol + assert abs(min(mesh.points[:, 0]) - 1.0) < tol + assert abs(max(mesh.points[:, 1]) - 2.0) < tol + assert abs(min(mesh.points[:, 1]) + 0.0) < tol + assert abs(max(mesh.points[:, 2]) - 3.0) < tol + assert abs(min(mesh.points[:, 2]) + 0.0) < tol + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 6.0) < tol + return + + +# # segfaults on travis, works locally +# def test_off(): +# pygalmesh.generate_from_off( +# 'elephant.off', +# 'out.mesh', +# facet_angle=25.0, +# facet_size=0.15, +# facet_distance=0.008, +# cell_radius_edge_ratio=3.0, +# verbose=False +# ) +# +# vertices, cells, _, _, _ = meshio.read('out.mesh') +# +# tol = 1.0e-3 +# assert abs(max(vertices[:, 0]) - 0.357612477657) < tol +# assert abs(min(vertices[:, 0]) + 0.358747130015) < tol +# assert abs(max(vertices[:, 1]) - 0.496137874959) < tol +# assert abs(min(vertices[:, 1]) + 0.495301051456) < tol +# assert abs(max(vertices[:, 2]) - 0.298780230629) < tol +# assert abs(min(vertices[:, 2]) + 0.300472866512) < tol +# +# vol = sum(compute_volumes(vertices, cells['tetra'])) +# assert abs(vol - 0.044164693065) < tol +# return + + +def test_extrude(): + p = pygalmesh.Polygon2D([[-0.5, -0.3], [0.5, -0.3], [0.0, 0.5]]) + domain = pygalmesh.Extrude(p, [0.0, 0.3, 1.0]) + pygalmesh.generate_mesh( + domain, "out.mesh", cell_size=0.1, edge_size=0.1, verbose=False + ) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-3 + assert abs(max(mesh.points[:, 0]) - 0.5) < tol + assert abs(min(mesh.points[:, 0]) + 0.5) < tol + assert abs(max(mesh.points[:, 1]) - 0.8) < tol + assert abs(min(mesh.points[:, 1]) + 0.3) < tol + assert abs(max(mesh.points[:, 2]) - 1.0) < tol + assert abs(min(mesh.points[:, 2]) + 0.0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 0.4) < tol + return + + +def test_extrude_rotate(): + p = pygalmesh.Polygon2D([[-0.5, -0.3], [0.5, -0.3], [0.0, 0.5]]) + edge_size = 0.1 + domain = pygalmesh.Extrude(p, [0.0, 0.0, 1.0], 0.5 * 3.14159265359, edge_size) + pygalmesh.generate_mesh( + domain, "out.mesh", cell_size=0.1, edge_size=edge_size, verbose=False + ) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-3 + assert abs(max(mesh.points[:, 0]) - 0.583012701892) < tol + assert abs(min(mesh.points[:, 0]) + 0.5) < tol + assert abs(max(mesh.points[:, 1]) - 0.5) < tol + assert abs(min(mesh.points[:, 1]) + 0.583012701892) < tol + assert abs(max(mesh.points[:, 2]) - 1.0) < tol + assert abs(min(mesh.points[:, 2]) + 0.0) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 0.4) < 0.05 + return + + +def test_ring_extrude(): + p = pygalmesh.Polygon2D([[0.5, -0.3], [1.5, -0.3], [1.0, 0.5]]) + edge_size = 0.1 + domain = pygalmesh.RingExtrude(p, edge_size) + pygalmesh.generate_mesh( + domain, "out.mesh", cell_size=0.1, edge_size=edge_size, verbose=False + ) + + mesh = meshio.read("out.mesh") + + tol = 1.0e-3 + assert abs(max(mesh.points[:, 0]) - 1.5) < tol + assert abs(min(mesh.points[:, 0]) + 1.5) < tol + assert abs(max(mesh.points[:, 1]) - 1.5) < tol + assert abs(min(mesh.points[:, 1]) + 1.5) < tol + assert abs(max(mesh.points[:, 2]) - 0.5) < tol + assert abs(min(mesh.points[:, 2]) + 0.3) < tol + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 2 * numpy.pi * 0.4) < 0.05 + return + + +# def test_heart(): +# class Heart(pygalmesh.DomainBase): +# def __init__(self, edge_size): +# super(Heart, self).__init__() +# return +# +# def eval(self, x): +# return (x[0]**2 + 9.0/4.0 * x[1]**2 + x[2]**2 - 1)**3 \ +# - x[0]**2 * x[2]**3 - 9.0/80.0 * x[1]**2 * x[2]**3 +# +# def get_bounding_sphere_squared_radius(self): +# return 10.0 +# +# edge_size = 0.1 +# d = Heart(edge_size) +# +# pygalmesh.generate_mesh( +# d, +# 'out.mesh', +# cell_size=0.1, +# edge_size=edge_size, +# # odt=True, +# # lloyd=True, +# # verbose=True +# ) +# +# return + + +def test_sphere(): + radius = 1.0 + s = pygalmesh.Ball([0.0, 0.0, 0.0], radius) + pygalmesh.generate_surface_mesh( + s, + "out.off", + angle_bound=30, + radius_bound=0.1, + distance_bound=0.1, + verbose=False, + ) + + mesh = meshio.read("out.off") + + tol = 1.0e-2 + assert abs(max(mesh.points[:, 0]) - radius) < tol + assert abs(min(mesh.points[:, 0]) + radius) < tol + assert abs(max(mesh.points[:, 1]) - radius) < tol + assert abs(min(mesh.points[:, 1]) + radius) < tol + assert abs(max(mesh.points[:, 2]) - radius) < tol + assert abs(min(mesh.points[:, 2]) + radius) < tol + + areas = compute_triangle_areas(mesh.points, mesh.cells["triangle"]) + surface_area = sum(areas) + assert abs(surface_area - 4 * numpy.pi * radius ** 2) < 0.1 + return + + +def test_halfspace(): + c = pygalmesh.Cuboid([0, 0, 0], [1, 1, 1]) + s = pygalmesh.HalfSpace([1.0, 2.0, 3.0], 1.0, 2.0) + u = pygalmesh.Intersection([c, s]) + + pygalmesh.generate_mesh(u, "out.mesh", cell_size=0.2, edge_size=0.2, verbose=False) + + mesh = meshio.read("out.mesh") + + vol = sum(compute_volumes(mesh.points, mesh.cells["tetra"])) + assert abs(vol - 1 / 750) < 1.0e-3 + return + + +if __name__ == "__main__": + test_sphere() -- 2.30.2